# Read in sheet 2 of our data
mtb_data <- read_excel(here::here('Data/mtb_stats.xlsx'), 'Sheet1')
mtb_data <- mtb_data %>% 
  # Clean up the label column
  mutate(label = str_replace_all(str_to_lower(label), '[:punct:]', ''),
         # Create a feature for the long-version of the names
         bike_category = case_when(
          label == 'tr' ~ 'Trail',
          label == 'xc' ~ 'Cross Country',
          label == 'dc' ~ 'Downcountry',
          label == 'am' ~ 'All Mountain',
          label == 'en' ~ 'Enduro',
          TRUE ~ 'red'
        ))

# Pull out the class labels
labels <- mtb_data %>% 
  select(label)


# Let's view the mtb_data output
# In any kable outputs, display NAs as blanks
opts <- options(knitr.kable.NA = "")

mtb_data %>% 
  drop_na('url') %>% 
  head(25) %>%
  # Fix up the headers by replacing the underscores with spaces
  rename_all(funs(str_replace_all(., "_", " "))) %>% 
  # Make everything proper capitalization
  rename_all(funs(str_to_title)) %>%
  kable() %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = F,
                font_size = 12) %>%
  # Make the header row bold and black so it's easier to read
  row_spec(0, bold = T, color = "black") %>% 
  scroll_box(height = "400px", width = "100%")

mtb_data %>% 
  group_by(bike_category) %>% 
  tally() %>% 
  arrange(desc(n)) %>% 
  # Start our visualization, creating our groups by party affiliation
  ggplot(aes(y = forcats::fct_reorder(bike_category, n), x = n)) +
  geom_col(fill = "slateblue", na.rm = T) +
  # Add a label by recreating our data build from earlier
  geom_label(aes(label = n),
             size = 5,
             # Scooch the labels over a smidge
             hjust = .25) +
  # Let's change the names of the axes and title
  xlab("Number of Bikes") +
  ylab("Category (label)") +
  labs(title = "Number of Mountain Bikes per Category")

mtb_data %>% 
  select(-label, -bike_category) %>% 
  DataExplorer::plot_bar(ggtheme = theme_classic(),
                         title = 'Distribution of Categorical Variables',
                         theme_config = theme(plot.title = element_text(hjust = 0, 
                                                                            color = "slateblue4", 
                                                                            size = 24),
                                                  plot.subtitle = element_text(hjust = 0, color = "slateblue2", size = 10),
                                                  plot.caption = element_text(color = "dark gray", size = 10, face = "italic"),
                                                  axis.title.x = element_text(size = 14),
                                                  axis.title.y = element_text(size = 14)),
                         maxcat = 15,
                         ncol = 2)

mtb_data %>% 
  select(-label) %>% 
  DataExplorer::plot_bar(ggtheme = theme_classic(),
                         by = 'bike_category',
                         by_position = 'fill',
                         title = 'Distribution of Categorical Variables',
                         theme_config = theme(plot.title = element_text(hjust = 0, 
                                                                            color = "slateblue4", 
                                                                            size = 24),
                                                  plot.subtitle = element_text(hjust = 0, color = "slateblue2", size = 10),
                                                  plot.caption = element_text(color = "dark gray", size = 10, face = "italic"),
                                                  axis.title.x = element_text(size = 14),
                                                  axis.title.y = element_text(size = 14)),
                         maxcat = 15,
                         ncol = 2)

DataExplorer::plot_density(mtb_data,
                             ggtheme = theme_classic(),
                             title = 'Distribution of Continuous Variables',
                             geom_density_args = list(fill = 'slateblue'),
                             theme_config = theme(plot.title = element_text(hjust = 0, 
                                                                                color = "slateblue4", 
                                                                                size = 24),
                                                      plot.subtitle = element_text(hjust = 0, color = "slateblue2", size = 10),
                                                      plot.caption = element_text(color = "dark gray", size = 10, face = "italic"),
                                                      axis.title.x = element_text(size = 14),
                                                      axis.title.y = element_text(size = 14)),
                             ncol = 3)

# Split data based on setting vs. no setting
no_setting <- mtb_data %>% 
  filter(is.na(setting))
setting <- mtb_data %>% 
  filter(!is.na(setting))



setting <- cbind(setting$model, setting$label, select_if(setting, is.numeric))
setting$model <- setting$`setting$model`
setting <- setting %>% select(-`setting$model`)
setting$label <- setting$`setting$label`
setting <- setting %>% select(-`setting$label`)

mean_by_setting <- aggregate(x=select(setting, -c(model, label)),
                             by=list(setting$model, setting$label),
                             FUN=mean)
mean_by_setting$model <- mean_by_setting$Group.1
mean_by_setting$label <- mean_by_setting$Group.2
mean_by_setting <- mean_by_setting %>% select(-c(Group.1, Group.2))

no_setting <- cbind(no_setting$model, no_setting$label, select_if(no_setting, is.numeric))
no_setting$model <- no_setting$`no_setting$model`
no_setting <- no_setting %>% select(-`no_setting$model`)
no_setting$label <- no_setting$`no_setting$label`
no_setting <- no_setting %>% select(-`no_setting$label`)

new_mtb_data <- data.frame(rbind(mean_by_setting, no_setting))

rownames(new_mtb_data) <- new_mtb_data$model

rm(no_setting)
rm(mean_by_setting)

mtb_correlation <- mtb_data %>% 
  # Get rid of price for now
  select(-price) %>% 
  # Select our variables of interest
  select_if(is.numeric) %>% 
  # Build our correlation matrix, such that missing values are handled by casewise deletion
  cor(use = 'complete.obs') 

# Convert our results into a tibble for easier manipulation
mtb_correlation_df <- mtb_correlation %>% 
  as_tibble() %>% 
  mutate(variable = colnames(mtb_correlation)) %>% 
  relocate(variable, everything())

# Build our correlation plot, using the angular order of the eigenvectors
corrplot(mtb_correlation,
         diag = F,
         col = COL2('PRGn'),
         tl.col = 'slateblue4',
         type = 'lower',
         method = 'color',
         order = 'AOE',
         title = 'Mountain Bike Feature Correlation'
         )

mtb_correlation_df %>% 
  pivot_longer(-variable, 
               names_to = 'correlated_variable', 
               values_to = 'correlation') %>% 
  filter(variable != correlated_variable) %>% 
  # Sort by the absolute value of correlation
  arrange(desc(abs(correlation))) %>% 
  filter((correlation > .90) | (correlation < -.90)) %>% 
  # Get rid of duplicative rows
  dplyr::distinct(correlation, .keep_all = T) %>% 
  pander()

# Impute missing values with column mean (not really best practice, but good enough)
for (c in 1:ncol(new_mtb_data)){
  if (is.numeric(unlist(new_mtb_data[,c]))){
    # print(colnames(new_mtb_data)[c])
    new_mtb_data[is.na(new_mtb_data[,c]), c] <- mean(unlist(new_mtb_data[,c]), na.rm=TRUE)  
  }
}

mtb_no_null <- new_mtb_data %>% 
                select(-price) %>%
                select_if(is.numeric) %>% 
                bind_cols(label = new_mtb_data$label) %>%
                drop_na()

mtb_pca <- prcomp(mtb_no_null %>% select(-label),
                  center = TRUE,
                  scale. = TRUE)

# Put our summary results into a dataframe - Justin switching this to cbind() works for me, not sure why
mtb_pca_df <- tibble(variable = c('Standard Deviation', 'Proportion of Variance', 'Cumulative Proportion')) %>% 
  cbind(summary(mtb_pca)$importance)


mtb_pca_df %>% 
  # Only display the first 6 columns
  select(c(PC1:PC5)) %>% 
  pander()

mtb_pca_df %>% 
  # Pivot our data so it's easier to visualize
  pivot_longer(-variable, 
               names_to = 'PC',
               names_prefix = 'PC') %>% 
  # Make the principal component column an integer so ggplot orders it from 1:17 properly
  mutate(PC = as.integer(PC),
         # Convert value to % (multiply by 100) so it's not a decimal
         value = 100*value) %>% 
  filter(variable == 'Proportion of Variance') %>% 
  ggplot(aes(x = PC, y = value)) +
  geom_point(size = 4, color = 'slateblue') +
  geom_line(alpha = .6, lwd = 2, color = 'slateblue') + 
  labs(title = 'Proportion of Variance Explained by Principal Components',
       x = 'Principal Component',
       y = 'Proportion of Variance (%)')

p_load(devtools,
       ggbiplot)

# jpeg('../Images/pca.jpg')

ggbiplot(mtb_pca,
              obs.scale = 1,
              var.scale = 1,
              groups = mtb_no_null$label,
              ellipse = TRUE,
              circle = FALSE,
              ellipse.prob = .5) + 
  theme(legend.direction = 'horizontal',
               legend.position = 'top')

# How many clusters are necessary? 4?

mtb_numeric <- mtb_no_null %>% 
  select(-label)
mtb_standard_scaled <- scale(mtb_numeric)

mtb_numeric <- mtb_no_null %>% 
  select(-label)

mtb_numeric <- mtb_no_null %>% 
  select(-label)

clusters <- 1:10
dists <- c()
for (c in 1:10){
  km <- kmeans(mtb_standard_scaled, centers=c, iter.max=1000)
  dists <- c(dists, km$tot.withinss)
}

# jpeg('../Images/Kmeans.jpg')
# plot(clusters, dists, type='l', xlab='Clusters', ylab='Total Sum of Squared Euclidean Distances')

# Plot our results
tibble(clusters = clusters,
       dists = dists) %>% 
  ggplot(aes(x = clusters, y = dists)) + 
  geom_point(size = 4, color = 'slateblue') +
  geom_line(alpha = .6, lwd = 2, color = 'slateblue') + 
  labs(title = "K-Means Clustering of MTB Data",
       subtitle = 'Method uses `tot.withinss` parameter to measure distances.',
       x = 'Clusters',
       y = 'Total Sum of Squared Euclidean Distances')

# Let's see where these clusters would end up on the 2D PCA plot
mtb_pca_scaled <- prcomp(mtb_standard_scaled,
                  center = F,
                  scale. = F)

pca_2_scaled <- as.matrix(mtb_standard_scaled) %*% as.matrix(mtb_pca_scaled$rotation[,1:2])

pca_km_scaled <- kmeans(pca_2_scaled, centers=3, iter.max=1000)


# Bring our PCA and k-means clusters results into our dataset
new_mtb_data %>% 
  cbind(pca_2_scaled) %>% 
  mutate(# Create a feature for the long-version of the names
         bike_category = case_when(
          label == 'tr' ~ 'Trail',
          label == 'xc' ~ 'Cross Country',
          label == 'dc' ~ 'Downcountry',
          label == 'am' ~ 'All Mountain',
          label == 'en' ~ 'Enduro',
          TRUE ~ 'red'
        ),
        # Bring our clusters in as a factor
        cluster = as.factor(pca_km_scaled$cluster)) %>% 
  # GG-plot our shit - lol
  ggplot() +
    geom_point(aes(x = PC1, 
                   y = PC2, 
                   color = cluster, 
                   shape = bike_category), 
               alpha = .9, 
               size = 3) +
    # Add our cluster centers in as well
    geom_point(data = as_tibble(pca_km_scaled$centers) %>%
                 mutate(cluster = as.factor(c(1, 2, 3))), 
               aes(x = PC1, 
                   y = PC2,
                   color = cluster), 
               shape = 10, 
               size = 7) + 
  # Color clusters accordingly
    scale_color_manual(values = c('slateblue4', 'gray', 'slateblue1'), name = 'Cluster') +
    labs(title = "K-Means Clustering of MTB Principal Components",
       subtitle = 'Assigned clusters denoted by color;\nBike categories denoted by shape;\nCluster centers denoted by large cross-hairs shape.',
       x = 'Principal Component 1',
       y = 'Principal Component 2')

#TODO let's look at this bottom cluster - both Niner bikes
# Niner has low reach numbers on its bikes - could be because we used the Medium for these!
# Based on PCA mapping, the blur tr, expic, Exie, Ripley, and Element all have less chainstay length, and less pistons?? wow, should we exclude piston count?? with more 2 piston bikes getting added, it evens out the average, so these aren't showing up as much anymore

p_load(ClusterR)


opt_gmm <- Optimal_Clusters_GMM(mtb_standard_scaled, 
                     max_clusters = 10, 
                     criterion = "BIC", 
                     dist_mode = "eucl_dist", 
                     seed_mode = "random_subset",
                     km_iter = 10, 
                     em_iter = 10, 
                     var_floor = 1e-10, 
                     verbose = T,
                     plot_data = T)

# Build our GMM model
mtb_gmm <- GMM(mtb_standard_scaled,
               gaussian_comps = 3,
               dist_mode = 'maha_dist', # Distance metric to use during seeding of initial means clustering
               seed_mode = 'random_subset', # How initial means are seeded prior to EM alg
               km_iter = 10, # Num of iterations of K-Means alg
               em_iter = 10, # Num of iterations of EM alg
               verbose = F
               )

# Run our predictions and convert it to a dataframe
mtb_gmm_pred <- predict_GMM(mtb_standard_scaled, 
                            mtb_gmm$centroids, 
                            mtb_gmm$covariance_matrices, 
                            mtb_gmm$weights) 

mtb_gmm_pred <- bind_cols(data.frame(mtb_gmm_pred$log_likelihood),
                       data.frame(mtb_gmm_pred$cluster_proba),
                       mtb_gmm_pred$cluster_labels) %>% 
  as_tibble()

# Rename our columns
names(mtb_gmm_pred) <- c('log_likelihood_c1', 
                         'log_likelihood_c2', 
                         'log_likelihood_c3',
                         'cluster1_prob',
                         'cluster2_prob',
                         'cluster3_prob',
                         'cluster_labels')

mtb_gmm_pred %>% 
  select(cluster_labels) %>% 
  bind_cols(mtb_no_null$label) %>% 
  # xtabs(~ ., data = .) %>% 
  table() %>% 
  bind_cols(cluster_labels = c('1', '2', '3')) %>% 
  relocate(cluster_labels, everything()) %>% 
  pander()

# Plot the densities for our clusters
cluster1_pred <-  mtb_gmm_pred %>% 
  filter(cluster_labels == 1) %>%
  ggplot(aes(x = cluster1_prob, y=..scaled..)) +
  # Plot the density and adjust the curvature so it looks better
  geom_density(position = 'stack', alpha = .1, fill = 'slateblue', color = 'slateblue', adjust = 10) +
  labs(title = 'Cluster 1',
       x = 'Probability',
       y = 'Density')

  
cluster2_pred <-  mtb_gmm_pred %>% 
  filter(cluster_labels == 2) %>%
  ggplot(aes(x = cluster2_prob, y=..scaled..)) +
  # Plot the density and adjust the curvature so it looks better
  geom_density(position = 'stack', alpha = .1, fill = 'slateblue', color = 'slateblue', adjust = 25) +
  labs(title = 'Cluster 2',
       x = 'Probability',
       y = 'Density')

cluster3_pred <- mtb_gmm_pred %>% 
  filter(cluster_labels == 3) %>%
  ggplot(aes(x = cluster3_prob, y=..scaled..)) +
  # Plot the density and adjust the curvature so it looks better
  geom_density(position = 'stack', alpha = .1, fill = 'slateblue', color = 'slateblue', adjust = 250) +
  labs(title = 'Cluster 3',
       x = 'Probability',
       y = 'Density')

gridExtra::grid.arrange(cluster1_pred, cluster2_pred, cluster3_pred, ncol = 3)

p_load(e1071,
       caret)

#convert all mountain category to enduro, dc -> Xc?
remap <- function(x, num){
  if (x=='am' || x=='en'){
    if (num){
      return(4)
    }
    else{
      return('Enduro')
    }
  }
  else if(x=='xc'){
    if(num){
      return(1)
    }
    else{
      return('Cross Country')
    }
  }
  else if(x=='dc'){
    if(num){
      return(2)
    }
    else{
      return('Downcountry')
    }
  }
  else if(x=='tr'){
    if(num){
      return(3)
    }
    else{
      return('Trail')
    }
  }
}
labels <- as.factor(unlist(lapply(new_mtb_data$label, remap, F)))


trainSVM <- function(x, y, idx, k='radial basis', g=0){
  xtest <- x[idx,]
  xtrain <- x[-idx,]
  ytest <- y[idx]
  ytrain <- y[-idx]
  
  if(k=='linearl'){
    clf <- svm(x=xtrain, y=ytrain, kernel=k)
  }
  else{
    if (g==0){
      clf <- svm(x=xtrain, y=ytrain, kernel=k)
    }
    else{
      clf <- svm(x=xtrain, y=ytrain, kernel=k, gamma=g)
    }
  }
  
  preds <- predict(clf, xtest)
  
  acc <- 0
  cm <- table(ytest, preds)
  for (i in 1:length(unique(labels))){
    acc <- acc + cm[i,i]
  }
  return(acc/sum(cm))
}



# Roughly 66.6% accuracy when treating down country as separate category
# But - roughly 68.8% accuracy when treating down country as XC, only 69% accuracy when treating downcountry as trail, suggests that downcountry bikes are more akin to trail than they are XC

folds <- createFolds(labels, k=10)


#Grid Search for SVM
kernels <- c('linear', 'polynomial', 'radial', 'sigmoid')
gammas <- seq(-5, 3, length.out=100)
gammas <- 10^gammas


#Change these below
X = pca_2_scaled
# X = mtb_standard_scaled
y = labels

results <- matrix(ncol=3, nrow=0)
colnames(results) <- c('acc', 'kernel', 'gamma')

for (k in kernels){
  if (k=='linear'){
    
    folds <- createFolds(labels, k=10)
    accs <- c()
    for (fold in folds){
      # print('here')
      acc <- trainSVM(X, y, fold, k='linear', g=0)
      accs <- c(accs, acc)
    }
    results <- rbind(results, c(mean(accs), k, 0))
  }
  else{
    for (g in gammas){
      # print('choosing gamma')
      # print(k)
      folds <- createFolds(labels, k=10)
      accs <- c()
      for (fold in folds){
        acc <- trainSVM(pca_2_scaled, labels, fold, k=k, g=0)
        accs <- c(accs, acc)
      }
    results <- rbind(results, c(mean(accs), k, g))
    }
  }
}

# results <- data.frame(results)
# results$acc <- as.numeric(results$acc)
# results$gamma <- as.numeric(results$gamma)
# results[which.max(results$acc), ]
# 
# 
# nonlinear_svm <- results[-(results$kernel=='linear'), ]
# 
# nonlinear_svm[which.max(nonlinear_svm$acc), ]
# On two axis, the best svm model is linear

## This cell only to visualize linear kernel
pcsvm <- svm(x=pca_2_scaled, y=labels, kernel='linear', gamma=0.2782559 )

dat <- data.frame(pca_2_scaled)

grid <- expand.grid(seq(min(dat[, 1]),max(dat[,1]),length.out=100),seq(min(dat[,2]),max(dat[,2]),length.out=100)) 
names(grid) <- names(dat)[1:2]
preds <- predict(pcsvm, grid)
df <- data.frame(grid, preds)

# ggplot(df, aes(x = PC1, y = PC2)) + 
#   geom_tile(aes(fill=preds)) +
#   geom_point(data = dat, aes(shape = labels), size = 2) + 
#   labs(title = "Support Vector Machine Classification",
#          x = 'Principal Component 1',
#          y = 'Principal Component 2')

ggplot(df, aes(x = PC1, y = PC2)) + 
  geom_tile(aes(fill = preds), alpha = .2) +
    scale_fill_manual(values = c('slateblue4', 'gray', 'slateblue1', 'green', 'red')) +
  geom_point(data = dat, aes(shape = labels, color = labels), size = 2) + 
  scale_color_manual(values = c('slateblue4', 'gray', 'slateblue1', 'green', 'red')) +
  labs(title = "Support Vector Machine Classification",
         x = 'Principal Component 1',
         y = 'Principal Component 2')

# DC as XC, using Linear
remap2 <- function(x, num){
  if (x=='am' || x=='en'){
    if (num){
      return(4)
    }
    else{
      return('Enduro')
    }
  }
  else if(x=='xc' || x=='dc'){
    if(num){
      return(1)
    }
    else{
      return('Cross Country')
    }
  }
  else if(x=='tr'){
    if(num){
      return(3)
    }
    else{
      return('Trail')
    }
  }
}
labels2 <- as.factor(unlist(lapply(new_mtb_data$label, remap2, F)))

pcsvm2 <- svm(x=pca_2_scaled, y=labels2, kernel='linear')

preds2 <- predict(pcsvm2, grid)
df2 <- data.frame(grid, preds2)

ggplot(df2, aes(x = PC1, y = PC2)) + 
  geom_tile(aes(fill=preds2)) +
  geom_point(data = dat, aes(shape = labels2), size = 2) + 
  labs(title = "Support Vector Machine Classification",
         x = 'Principal Component 1',
         y = 'Principal Component 2')